function [dof_map, V, T, u1, u2, err] = Cook_Membrane(n_round, Young, nu, FE_Order, method, FE_Type)
% function [dof_map, V,T, u1, u2, err] = Cook_Membrane(n_round, Young, nu, FE_Order, method, FE_Type)
%%
%
% benchmark problem: Cook's membrane example, see
%      K.S. Chavan et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007)
%      4075-4086.
% and its ref
%  [1]
%  [10]
%  [27]
%
%  Test pass on Aug, 2012
%
%  Author: Dr. Xian-Liang Hu
%

if nargin < 2
    % default physical parameters for elastic problems
    Young = 250;   % N/mm^2
    nu = 0.3;     % limit is 0.5
end

if nargin < 4
    FE_Order = 3;   % default polynomial order.
end

if nargin < 5
    method = 2;  % using FEM style by default
end

if nargin < 6
    FE_Type = 'BB';  % default FE basis using Bernstein Polynomial
end

Quad_Order = 13;

%%%%%%%%%%%%
% generate initial mesh
% tt = [0;0.5;1];
% [V,T] = init_mesh_rectangle(tt, tt, 'type3A');
% [T, E, ET, TE] = build_fem_mesh(V, T);
% for round= 1:n_round
%     [V,T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET);
% end

load MembraneMesh8500 V T E ET TE;

% modify the corner, so as the domain
x = V(:,1); y = V(:,2);
VA = [0, 0]; VB = [48, 44]; VC = [48, 60]; VD = [0, 44];
V = (1 - x).*(1-y)*VA + x.*(1-y)*VB + x.*y*VC + (1-x).*y*VD;

% ZERO=zeros(size(V,1),1); trisurf(T,V(:,1),V(:,2),ZERO);view(0,90);
fprintf('There are %d triangles.\n', size(T,1));


%%%%%%%%%%
%  
bdr_Dirichlet = find((V(E(:,1),1) < 1e-5) & (V(E(:,2),1) < 1e-5));
bdr_Neumann = find((V(E(:,1),1) > 48 - 1e-5) & (V(E(:,2),1) > 48 - 1e-5) & ...
         (V(E(:,1),2) < 53) & (V(E(:,1),2) > 51));

[dof_map, u1, u2, t] = elastic_fem(V, T, TE, ET, FE_Type, FE_Order, Quad_Order, method, ...
                                     Young, nu, @fun_f, bdr_Dirichlet, @fun_ub, bdr_Neumann, @fun_surface_f);
           
%%%%%%%%%%%%%%%%%%%%%%
%%% this is for the linear FE case:
% n_dof = max(max(dof_map)); 
% u = zeros(n_dof,1); u(dof_map) = u1;
% trisurf(T,V(:,1),V(:,2),u);
% n_dof =max(max(dof_map));
% u = zeros(n_dof,1); u(dof_map) = u1;
% v = zeros(n_dof,1); v(dof_map) = u2;
% ONES = ones(size(V,1),1);
% fac = 0; trisurf(T,V(:,1) + fac*u,V(:,2) + fac*v,fac*ONES.*V(:,1));view(0,90); % original plate
% hold on;
% fac = 0.5; trisurf(T,V(:,1) + fac*u,V(:,2) + fac*v, fac*ONES.*V(:,1));view(0,90);% plus half of the displacement
% fac = 1; trisurf(T,V(:,1) + fac*u,V(:,2) + fac*v, fac*ONES.*V(:,1));view(0,90); % plus the displacement
% axis([0 50 0 100]);
% hold off;


%%%%%%%%%%%%%%%%%%%%%%
% there is no analytic solutions in this case
err = 1;

%%%%%%%%%%%%%%%%%%%%%%
%% order d finite element solutions:
disp_d = FE_Order; tri_temp = template_mesh_tri(disp_d);
[u_h, Tris1, Points1] = fe_solution_bb(V,T, u1, FE_Order, tri_temp, disp_d);
[v_h, Tris2, Points2] = fe_solution_bb(V,T, u2, FE_Order, tri_temp, disp_d);


%%%%%%%%%%%%%%%%%%%%%%
% % visualization and print information
ONES = ones(size(Points1,1),1);
fac = 0; trisurf(Tris1,Points1(:,1) + fac*u_h, Points1(:,2) + fac*v_h, fac*ONES.*Points1(:,1));view(0,90); % original plate
hold on;
% fac = 0.5; trisurf(Tris1,Points1(:,1) + fac*u_h, Points1(:,2) + fac*v_h, fac*ONES.*Points1(:,1));view(0,90);% plus half of the displacement
fac = 1; trisurf(Tris2,Points2(:,1) + fac*u_h, Points2(:,2) + fac*v_h, fac*ONES.*Points2(:,1));view(0,90); % plus the displacement
axis([0 50 0 100]);
hold off;


% subplot(1,2,1); trisurf(Tris1, Points1(:,1), Points1(:,2), u_h);  %plot error for u
% subplot(1,2,2); trisurf(Tris2, Points2(:,1), Points2(:,2), v_h);  %plot error for v


% fprintf('The inf norm or error(u) = %e,  error(v) = %e.\n',  erru, errv);

end


%%%%%%%%%%%%%%%%%%%%%%%%
%  this is boundary condition
% which is clamped on the left bounary
% 
function [u_1, u_2] = fun_ub(xx, yy, varargin)
    u_1 = 0*xx;
    u_2 = 0*yy;
end

%%%%%%%%%%%%%%%%%%%%%%%%
%  the right hand side of the elastic problems, means the force on the
%  elements
% 
function [u_1, u_2] = fun_f(xx, yy, varargin)
    u_1 = 0*xx;
    u_2 = 0*yy;
end


%%%%%%%%%%%%%%%%%%%%%%%%
%  The surface force applied on boundary, we referred it as Neumann boundary
%  which is only loaded on the right boundary 
%  with 100 N/m in positive y direction
%
function [f_1, f_2] = fun_surface_f(xx, yy, varargin)
    f_1 = 0*xx;
    f_2 = 0*yy;
    idx = (xx > 48 - 1e-5);
    f_2(idx) = 100;
end